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We study s-wave superconductivity in the two dimensional attractive Hubbard model in an applied magnetic 
field, assume the extreme Pauli limit, and examine the role of spatial fluctuations in the coupling regime cor¬ 
responding to BCS-BEC crossover. We use a decomposition of the interaction in terms of an auxiliary pairing 
field, retain the static mode, and sample the pairing field via Monte Carlo. The method requires iterative solution 
of the Bogoliubov-de-Gennes (BdG) equations for amplitude and phase fluctuating configurations of the pairing 
field. We establish the full thermal phase diagram of this strong coupling problem, revealing T c scales an order 
of magnitude below the mean field estimate, highlight the spontaneous inhomogeneity in the field induced mag¬ 
netization, and discover a strong non monotonicity in the temperature dependence of the low energy density of 
states. We compare our results to the experimental phase diagram of the imbalanced Fermi gas at unitarity. This 
paper focuses on the magnetized but homogeneous (breached pair) superconducting state, a companion paper 
deals with the Fulde-Ferrell-Larkin-Ovchinnikov (FFFO) regime. 


I. INTRODUCTION 


For an electron system in a superconducting state the 
Meissner effect characterizes the response to a magnetic field. 
In type II superconductors there is flux penetration beyond a 
threshold h c \ in the form of an Abrikosov latticd 1 2 , before su¬ 
perconductivity (SC) is finally lost at the ‘orbital critical field’ 
h° 2 b . The magnetic field also couples to the spin of the elec¬ 
trons, and tends to break an ‘tl’ pair. This effect is detri¬ 
mental to SC, and, if orbital effects were irrelevant, SC order 
would be lost at some ‘Pauli limiting’ field, h^ 2 , say 3 A The 
ratio of these critical fields, a = h° 2 b /h^ 2 , defines the Maki 
parameter and is roughly® Ao /cp, where Ao is the zero tem¬ 
perature gap in the SC state and ep is the Fermi energy. 

In most superconductors a <C 1, so the Pauli suppression 
effects never show up. There are however three scenarios 
where it becomes relevant, (a) If ep is suppressed strongly 
by correlation effects, as in heavy fermions where the sup¬ 
pression factor can be ~ 10 3 due to Fermi liquid corrections 8 , 
(b) for two dimensional systems, the layered organics, say, or¬ 
bital effects are irrelevant for an ‘in plane’ field, and (c) for 
neutral Fermi gases, as in cold atomic systems, the magnetic 
effects would be related only to spin. Recent discoveries on 
the heavy fermion 9 17 CeCoIns, the ft-BEDT based layered 
superconductors 18 25 , iron pnictidesPi® and population im¬ 
balanced cold Fermi gases 31 40 , make the Pauli limit relevant. 

Early extension^ of the BCS scheme to finite Zeeman 
field (neglecting orbital effects) predicted that, in the con¬ 
tinuum, the superconducting T c decreases with applied field 
up to a critical value, hi, say, and the thermal transition re¬ 
mains second order. Beyond hi, one would have expected a 
SC state with a first order thermal transition, but the ground 
state actually becomes modulated, in the spirit predicted by 
Fulde and Ferrell (FF j^I and Larkin and Ovchinnikov (LOj^J 
and stays so till all order is lost at some h c 2 . The FFLO state 
is characterized by periodic spatial modulation of the super¬ 
conducting order parameter and magnetization. The modu¬ 
lations in these quantities are complementary and the nodes 
of the superconducting order correspond to the maximum of 


the magnetization, and vice versa. The zero temperature 
and h < hi system is an ‘unpolarised superfluid’ (USF), its 
finite temperature counterpart is a ‘breached pair’ (BP) state 
and the hi < h < h c 2 window is FFLO. The BP phase is the 
finite temperature extension of the USF, with homogeneous 
superconducting order spatially coexisting with finite uniform 
polarization. The original scenario does not support any first 
order transition between the BP phase and the normal state. 

Experiments bear out some features of this scenario with 
solid state studies focused on probing the FFLO state while 
cold atom experiments probe the general effect of imbalance 
on pairing. The FFLO signatures in CeCoIn 5 include specific 
heat 9 , magnetic torque 44 , muon_spin relaxation^, NMR 1116 , 
and magnetic neutron scattering 15 while in the ft-BEDT based 
organics there is indirect evidence®^ for a modulated state 



FIG. 1. Color online: (a) Comparison of T c scales obtained from 
the mean field calculation (upper curve) and our static auxiliary field 
(SAF) Monte Carlo technique (lower curve). In the SAF data BP-II 
represents a breached pair state that undergoes a second order tran¬ 
sition to the partially polarized Fermi liquid (PPFF), while BP-I un¬ 
dergoes a first order thermal transition to the PPFF. Beyond BP-I the 
system exhibits FFLO order up to some critical field, (b) Polarization 
-vs- temperature phase diagram inferred from the SAF calculation, 
plotted in the spirit of the experimentally obtained unitary Fermi gas 
result 32 . The ‘unstable’ region is phase separated. 










at large in-plane fields. In KFe 2 As 2 also thermal expansion 
and magnetostriction suggests the occurrence of Pauli limited 
superconductivit)®. For cold atoms, fermionic superfluidity 
with population imbalance has been probed in detail with the 
Fermi gas tuned to unitarit} 32 revealing an ‘universal’ phase 
diagram. 

The microscopic models for superconductivity (or superflu¬ 
idity) in these systems are widely different but they share the 
features of (i) a ‘homogeneous’ magnetized superfluid state 
near T c at intermediate fields, (ii) a possible FFLO state at 
higher fields, and (iii) being in a coupling regime well be¬ 
yond the reach of mean field theory (at least for the atomic 
gases). Taking (iii) as our point of departure we address these 
issues by studying the Zeeman field dependence in the at¬ 
tractive two dimensional Hubbard model at intermediate cou¬ 
pling, U/t = 4 (see later). This corresponds roughly to the 
maximum T c in the BCS-BEC crossover window, and cru¬ 
cially involves amplitude and phase fluctuations in describing 
the thermal physical 6 49 . Our main results, using a recently 
developed Monte Carlo (MC) approach, are the following: 

1. We observe that in the imbalanced problem, as in the 
case of balanced Fermi gaseJSKI, fluctuation ef¬ 
fects suppress T c scales by a factor of more than 4 com¬ 
pared to widely used mean field theory. 

2. Intermediate fields allow for a temperature window over 
which the superfluid supports significant magnetization 
which, although homogeneous on the average, shows 
noticeable configurational fluctuation. 

3. At high fields the superfluid shows a first order tran¬ 
sition to the normal state on heating, but cooling in 
this field window inevitably traps the system into a 
metastable FFLO state. 

4. The spin resolved density of states shows a pseudogap 
(PG) feature that is strongly non monotonic in temper¬ 
ature: the pseudogap weakens initially with increas¬ 
ing temperature and then deepens again beyond a scale 
Tmax • The applied field dramatically suppresses T ma;r . 

We characterize the thermal state via real space maps, 
the structure factors associated with the superfluid and mag¬ 
netic order, the spin resolved momentum distribution of the 
fermions, and the density of states. 

The rest of the paper is organized as follows. In Section- 
II we discuss the model and the methods used to study it. 
Section-Ill contains the results. Section-IV discusses possi¬ 
ble limitations of our numerical scheme, suggests the con¬ 
nection to Ginzburg-Landau phenomenology, and relates our 
predictions to some cold atom experiments. Section-V con¬ 
cludes with our key observations. An appendix describes the 
Hubbard-Stratonovich transformation and the related approx¬ 
imations in detail. 


II. MODEL AND METHOD 
A. Model 

We study the attractive two dimensional Hubbard model 
(A2DHM) on a square lattice in the presence of a Zeeman 
field: 

h = h 0 - - \u\ (!) 

i i 

with, H 0 = where Uj = -t only 

for nearest neighbor hopping and is zero otherwise, gi z = 
(1/2) (n^ — riif). We will set t = 1 as the reference energy 
scale, ti is the chemical potential and h is the applied mag¬ 
netic field in the z direction. U > 0 is the strength of on-site 
attraction. We will use U/t = 4. 

We wish to explore the physics beyond weak coupling, i.e, 
short coherence length. This requires retaining fluctuations 
well beyond mean field theory (MFT). We accomplish that 
as follows. We use a ‘single channel’ Hubbard-Stratonovich 
(HS) decomposition of the interaction term in terms of an aux¬ 
iliary complex scalar field A i{r) = \Ai(r)\e iei ^ T \ A com¬ 
plete treatment of the resulting problem handles the full spa¬ 
tial and imaginary time, (i,r), dependence of the A - this 
is possible only within quantum Monte Carlo - while mean 
field theory imposes a spatially periodic pattern and ignores 
the r dependence. We ignore the ‘time’ dependence of the 
A, but completely retain the spatial dependence. This, as we 
shall see, makes our method ‘mean field’ at zero temperature, 
T = 0, but retains the crucial thermal fluctuations of the am¬ 
plitude and phase of A i that control T c scales, etc. We discuss 
the formal structure of this approximation in detail in the Ap¬ 
pendix, and its limitations in the Discussion section. 

The static A^ problem is described by the coupled effective 
Hamiltonian: 

H eff = H 0 - h a iz + Yx*4+4i + h - c ) + H <* 

i i 

where H c i = JA is the stiffness cost associated with the 
now ‘classical’ auxiliary field. The equation above indicates 
how the fermions see the pairing field. The pairing field con¬ 
figurations in turn are controlled by the Boltzmann weight: 

P{Ai} <x Tr CjCt e-0 H «" (3) 

This is related to the free energy of the fermions in the con¬ 
figuration {A^}. For large and random A^ the trace needs to 
be computed numerically. We generate the equilibrium {A^} 
configurations by a Monte Carlo technique (see later) diago- 
nalising the fermion Hamiltonian H e ff for every attempted 
update of the auxiliary fields. 

B. Numerical method: Monte Carlo and variational 
calculation 

Mean field theory has been the standard tool for exploring 
the effect of a Zeeman field on the superconductor. However, 






even though MFT may be reasonable in capturing the ground 
state, inclusion of amplitude and phase fluctuations is essen¬ 
tial as one moves beyond the U/t <C 1 window. This issue 
has been widely discussecP® i n the context of the zero field 
BCS to BEC crossover. 

Fluctuation effects have been found to suppress the T c , 
compared to MFT estimates, both at intermediate and strong 
coupling. Measurements on the 3D unitary gas indicateJ 70 71 
a peak T C /E F ~ 0.167, while various theoretical estimates 
at unitarity include (a) a mean field based result 5 ^ suggest¬ 
ing T c /E f ~ 0.6 (b) a T-matrix based result 51 suggesting 
T C /E F ~ 0.16, (c) fluctuation corrected mean field theory22, 
yielding TJE F ~ 0.245, and (d) Monte Carlo estimates 
y ieldin^pKS T c / E F ~ 0.15 — 0.25. A very recent experiment 
on a 2D cold Fermi gas indicated 72 ^ a peak T C /E F ~ 0.16, 
while an interpolative theory estimate suggests T C /E F ~ 0.1. 
Results on the 2D Hubbard model indicate 73 T c /t = 0.16. 
Corrections beyond mean field theory, it is obvious, are es¬ 
sential for an accurate description beyond weak coupling. 

We include thermal fluctuations via our static auxiliary field 
(SAF) scheme, which, implemented using Monte Carlo, can 
access system sizes larger than typical quantum Monte Carlo 
(QMC) calculations. This has several advantages: (i) it pro¬ 
vides an accurate estimate of the T c , (ii) at high fields it helps 
in accessing spatially modulated (FFLO) paired states which 
may have a large wavelength, and (iii) it allows calculation of 
dynamical properties without the need for any analytic con¬ 
tinuation. 

In order to make the study numerically less expensive the 
Monte Carlo is implemented using a cluster approximation, 
in which instead of diagonalising the entire L x L lattice for 
each local update of the A^ a smaller cluster, of size L c x 
L c , surrounding the update site is diagonalised. We mostly 
used L = 24 and L c = 6 for the results in this paper. The 
cluster approximation has been extensively benchmarked, and 
used successfully in the zero field case 49 . We will discuss the 
limitations of the SAF approach and cluster based update at 
the end of the paper. 

At zero temperature within the SAF scheme the energy is 
minimized over static configurations of the field A F We have 
carried out variational calculations at several fixed values of y, 
at different ft, exploring the following kinds of periodic con¬ 
figuration: (i) ‘axial stripes’: A^ ~ A 0 cos(go^), and diago¬ 
nal stripes A * ^ Ao cos (q(xi Eyi)), and (ii) two dimensional 
modulations, A i ~ A 0 (cos (q.Xi) + cos (q.yi)), and of course 
(iii) the unpolarised superfluid (USF) state A^ = Ao. We min¬ 
imize the energy with respect to the q , and A 0 (assumed real). 
This paper focuses on the uniform state, the FFLO regime is 
discussed in detail in a companion paper. 


C. Parameter regime and indicators 

Any real space numerical calculation requires a system with 
linear dimension L ^ o, where £o is the T = 0 coherence 

length, to accurately capture the SC state. Since £o increases 
with reducing U/t, this puts a limit on the U/t window that 
we can explore. The results in this paper are at U = 41, both 


within Monte Carlo and the variational scheme. We have also 
explored U = 2t variationally but it requires L ~ 48 to access 
modulated phases so we have not been able to do MC in that 
regime. At U/t = 4 we have explored the h—T dependence at 
multiple values of y below half-filling (the physics above half¬ 
filling can be inferred from this) but the qualitative features are 
independent of the choice of y so this paper focuses on y m 
—0.2 1, where the density is n « 0.94 (independent, roughly, 
of ft or T). We have studied the temperature dependence at a 
large number of fields in the window h/t ~ [0 : 1.5]. Beyond 
the global features of the ft — T phase diagram, we will discuss 
three field values, typical of three response regimes. 

We use the following indicators to characterize the sys¬ 
tem: (i) Monte Carlo snapshots of |A^|, the phase correlation 
cos(0q — 0i) where is the angle at a fixed reference site 
on the lattice, the magnetization variable m* = (n^ — n 4 ), 
and particle number rii = (n^ + n^). These explicitly high¬ 
light the spatial fluctuation with increasing temperature, and 
the modulated nature in the FFLO window, (ii) We keep track 
of the structure factors, Sa (q) and S m ( q), defined as: 

i,3 

Sm^) = Ej2(m i m j )e i ^-^ 

i,3 

where, N = L 2 . (iii) We monitor the bulk magnetization 
and the SC order parameter, Sa (q = 0; T, ft). (iv) We com¬ 
pute the momentum occupation number ((n^ a )) that carries 
the signature of imbalance and FFLO modulation. Finally, 
(v) we compute the spin resolved and total fermionic density 
of states (DOS). 

III. RESULTS 

In what follows we first highlight the huge difference be¬ 
tween the mean field results and that of our Monte Carlo ap- 



FIG. 2. Color online: The variational estimate of energy for varying 
amplitude A 0 and different modulation vectors q = {2n7r/L,0}. 

(a) h = 0.50t, where the ground state is USF, i.e, q = (0,0) and 

(b) h = 0.95 1 where the ground state is amplitude modulated. We 
have shown data only for axial modulations, the actual comparison 
is made for the full set described in the text. 



















FIG. 3. Color online: Ground state [i — h phase diagram obtained 
from the variational scheme, showing the unpolarised superfluid 
(USF), modulated (LO) and partially polarized Fermi liquid (PPFL) 
regions. There is no homogeneous superfluid state with finite mag¬ 
netization, i.e, BP, at T — 0. 


proach due to that of thermal fluctuations in this coupling 
regime. We then take a step back to illustrate the working 
of the variational approach to the ground state and the /i — h 
phase diagram that emerges. Following this we move on to 
a detailed discussion of thermal properties, in particular the 
difference between ‘cooling’ and ‘heating’ the system, sug¬ 
gestive of the presence of metastable states. We show detailed 
results for what we feel are three broad field regimes: (i) Weak 
field, where the T c is only modestly modified with respect to 
h = 0, the thermal transition is second order, and there is 
hardly any magnetization for T <T C . (ii) Intermediate field, 
where T c is noticeably lower, the thermal transition is still sec¬ 
ond order, but there is a window ST = T c — T > 0 where the 
system simultaneously shows superfluid order and magneti¬ 
zation, characteristic of the ‘breached pair’ state, (iii) Strong 
field, where the SC shows a first order thermal transition, and 
there is a metastable FFLO state over a wide temperature win¬ 
dow. 

Fig.l presents the primary contrast between MFT and the 
MC result. Fig.l.(a) presents the h — T phase diagram indi¬ 
cating regions of first and second order thermal transition and 
the regions of BP and FFLO character. A much more detailed 
phase diagram will be shown in Fig.4. 

Fig.l.(b) shows the MC phase diagram in terms of the in¬ 
ferred magnetization and temperature to create a parallel with 
cold fermionic system^ where the physics is probed for 
a fixed population imbalance (“magnetization”) rather than 
a fixed applied field. At T = 0 the entire USF window, 
0 < h < 0.85 t, collapses to the origin, and the first order 
jump to the LO state involves a magnetization discontinuity 
m ~ 0.28. Magnetization in the LO ground state is up to 
m ~ 0.37 beyond which we have the ‘normal’ partially po¬ 
larized Fermi liquid. The BP state is the finite T extension 
of the USF and occupies a widening window as T increases 
and then shrinks again as T —>> T®. The ‘unstable’ region 
is the magnetization discontinuity between the high tempera¬ 
ture PPFL state and the low T nearly unpolarised state in the 


1st order transition window. The LO T c ’s are small and the 
LO phase occupies a small low temperature sliver in the large 
m region. This picture helps understand the cold atom experi¬ 
ments where the population imbalance, rather than the applied 
field, is the primary variable. We will take up this comparison 
at the end of the paper. Fluctuations suppresses the T c to well 
below the MF value, as observed earlier in balanced Fermi 
gases 51 55 . The presence of imbalance (or an applied field) 
suppresses the T c more rapidly. 

We have used the variational approach described earlier to 
determine the ground state, in the same spirit as Chiesa et 
al . 747 5, wherein diagonal, uniaxial and checkerboard patterns 
of A i were compared to determine the ground state. 

Fig.2.(a) and Fig.2.(b) shows the dependence of the energy 
on the ‘magnitude’ A 0 , of the pairing field, for several values 
of q. Panel (a) is for intermediate field, h = 0.5£, where the 
ground state is still homogeneous, i.e, at q = (0, 0). Panel (b), 
at h — 0.95t shows an absolute minimum at q = (7t/3, 0), an 
axial Larkin-0vchinnikov state. 

The variationally determined fi — h phase diagram is shown 
in Fig.3. At low h the system is a homogeneous unmagnetised 
superfluid (USF). One may have expectecf 3 this to undergo 
a transition to a partially polarized Fermi liquid (PPFL) at a 
field h c = Ao/y / (2), the naive Pauli limit. However, as pre¬ 
dicted by Fulde and Ferrell 42 and Larkin and Ovchinnikov^, 
and confirmed by several later studies, we find that a A^ mod¬ 
ulated state with finite magnetization intervenes between the 
USF and the PPFL. We designate the USF to LO transition as 
h c i and the LO to PPFL transition as h c 2 . Both these fields 
increase with /i. We will discuss the detailed behavior within 
the LO window elsewhere. 


A. Overview of thermal phase diagram 

Mean field theory for 5-wave superconductors in a mag¬ 
netic field indicate that (in the continuum case) the normal to 
SC thermal transition continues to be second order from h — 0 
to a finite field, beyond which the system shows a first order 
transition, but now to a modulated superfluid phase 42 43 . 

The simultaneity of the second to first order change and 
transition from the q = (0, 0) to a finite q state is probably 
specific to continuum mean field theory. Additionally, the MF 
prediction of T c scales, etc , is valid only in the weak coupling 
limit. 

In the presence of an underlying lattice, even MFT suggests 
a field window over which one can have a first order SC to nor¬ 
mal transition, see Fig.l, although the transition temperature 
is badly overestimated. Beyond another higher field the lattice 
based MFT predicts a modulated state. 

Fig.4 shows the phases revealed by heating from the mean 
field ground state (left) and cooling (right) from a disordered 
high temperature state. The thermal transition from the SC to 
normal state is second order up to a field hi ~ 0.7 1 beyond 
which it becomes first order (with the ordered state still being 
at q = (0, 0)). For h < 0.7£ the results are path independent 
but for 0.7f < h < 0.85 1 the system gets trapped in a LO 
state on cooling although the ground state is still USF. Beyond 
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FIG. 4. Color online: The field-temperature, h — T, behavior suggested by (a) heating from the variational ground state, and (b) cooling from 
a random high temperature state. We show the phases that emerge, the T c scales, as well as the dominant fluctuation in the disordered regime 
(following a convention described in the text). In both panels the change from a second order to first order BP to PPFL thermal transition 
occurs consistently at h ~ 0.7 1. In the heating process the first order BP-PPFL transition encounters a region with strong LO fluctuations. 
This regime shows LO fluctuations on cooling as well and the system remains trapped in a fragmented LO state, rather than transit to the BP 
phase, as T is lowered. The USF to LO transition in the ground state occurs at h ~ 0.85 1. The ‘LO’ window in (a) refers to the genuine 
ground state, while in (b) it also includes the metastable LO region. 


h ~ 0.85 1, where the ground state is LO the results are again 
path independent. 

In the first order transition window, hi < h < h c i, the 
USF ground state thermally evolves into BP at finite T and 



FIG. 5. Color online: (a)-(b) Monte Carlo results for the temperature 
dependence of (a) 5(0, 0) and (b) the magnetization, for heating 
(solid line) and cooling (open circle), (c)-(d). Mean field results on 
the order parameter (c), and the magnetization (d). Note that the T 
range in (c)-(d) is about 5 times larger than in (a)-(b). Also, the MC 
based order parameter shows an almost linear drop with T at low 
temperature while the MF order parameter is expectedly flat. The 
magnetization results, (b)-(d), despite their overall similarity differ 
in the low h,T>T c window. 


then shows a transition to a PPFL state where the fluctuations, 
surprisingly, have LO character. On cooling down from a dis¬ 
ordered state the system fails to attain a q = (0,0) state and 
instead shows strong LO signatures. This MC inferred LO 
state is energetically higher than the variational USF state so 
this is a sign of metastability. We would characterize this state 
in terms of the various indicators in a later section. 


1. Fluctuation regime 

While long range order is only observed for T < 0.2 1, 
we wanted to probe if there is a significant window above T c 
where fluctuation effects of q = (0,0) or finite q pairing can 
be seen. We define the cut off to the fluctuation regime as the 
temperature at which the ratio between the highest magnitude 
of the structure factor peak to that at the neighboring k-point 
is ~ 1.5. The regimes of strong BP fluctuation and strong LO 
fluctuation are marked in Fig.4 in this spirit. 


2. Thermodynamic properties 


Fig.5 shows the thermal evolution of q = (0,0) structure 
factor peak, 5(0,0), and magnetization m(T) for the mag¬ 
netic fields characteristic of the low, intermediate and high 
field regimes. The two panels on top show the MC based re¬ 
sults while the lower panels are based on MFT. The MC and 
MFT results have a gross similarity but (i) the MF T c scales 
are four times larger, (ii) even on a normalized, T/T c scale, 
the MC order parameter shows a quicker drop with temper¬ 
ature associated with the 0(2) nature of the superconduct¬ 
ing problem, while the MF plot is much flatter due to ab¬ 
sence of phase fluctuations, (iii) the MF magnetization has 
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FIG. 6. Color online: Thermal evolution of the superfluid (5a (q)) and magnetic (5 m (q)) structure factor at h — 0.2t. By the time the 
magnetization picks up a reasonable value (extreme right) the superfluid order has been lost. 


dm/dT < 0 above T c , while the MC result, up to intermedi¬ 
ate field, clearly shows dm/dT > 0. 

The MC data also reveal that in the high field region, cor¬ 
responding to a first order transition, the USF ground state 
is not recovered on cooling and the response becomes his¬ 
tory dependent. If one heats from the USF ground state to 
the T > T c and cools again to T = 0 a large magnetization 
state is obtained! This state, as we will see, has prominent LO 
correlations. 

In what follows we provide a detailed description of the 
thermal response of the imbalanced superconductor for three 
typical field regimes. 


B. Low field response: the unpolarised superfluid 

We characterize 0 < h < 0.3t roughly as the Tow field’ 
regime. The T c is still within 10% of the h = 0 value and 
the magnetization even near T c is quite small (ra(T c ) ^0.02 
at h = 0.3t). As representative of this regime we show data 
for h = 0.2t in Fig.6, where the upper panel is the pairing 
structure factor Sa (q) and the lower panel is the magnetic 
structure factor 5 m (q). Sa loses its ordering feature at T ~ 
0.16t where the 5 m (0, 0) still remains < 10 -5 . 

There are no finite q features in the magnetic structure fac¬ 
tor. The thermal transition is reversible and no thermal history 
effects show up. 


C. Intermediate field: breached pair state 

Next we consider the intermediate field regime of 0.3t < 
h < 0.7t. Over this window T c (h ) falls significantly and the 
magnetization below T c reaches ~ 0.15 (at h = 0.7t). Fig.7 
shows spatial features at h = 0.5t. While the ground state 
is still a homogeneous USF the increase in T leads quickly 
to development of finite magnetization. ‘Unpaired’ fermions 


coexist with a q = (0,0) condensate. This is a breached pair 
state. 

We characterize this phase in Fig.7 through the T depen¬ 
dence of the following indicators: (a) the pairing ampli¬ 
tude \A(x,y)\, (b) phase correlation cos(#o — 0 x , y ) where 
#o is the phase at a reference site, (c) pairing structure factor 
5 a (q), (d) magnetization m(x,y), (e) magnetic structure fac¬ 
tor 5 m (q), and (f) number density n(x, y ). (a), (b), (d) and (f) 
are for a single MC snapshot, while (c) and (e) are thermally 
averaged. The T c in this case is ^ 0.13t. 

With increase in temperature MC snapshots indicate that 
the \A(x,y)\ becomes inhomogeneous (although a thermal 
average would be homogeneous again), and the phases be¬ 
gin to decohere. The T dependence of |A|, phase correla¬ 
tions, and the SC structure factor are not qualitatively differ¬ 
ent from what we see at weak field but m(x,y) shows a depar¬ 
ture. Between T = 0.lit and 0.15t, i.e , across T c , we observe 
the emergence of significant magnetization in ‘clumps’. The 
magnetization, crudely, follows a pattern that is spatially com¬ 
plementary to the SC order. The local magnetization can reach 
a value ~ 0.4 even for T <T C (the system average however 
is much smaller). The 5th row shows the magnetic structure 
factor, essentially a diffuse peak around q = (0, 0), while the 
last row shows the density profile (almost homogeneous). 

We have calculated the momentum occupation number 
n a ( k) = «4 a c k(T ». In Fig.8 we show n t (k) and n l( k ) 
at h = 0.5t for different temperatures. At low temperature 
where the system is unpolarised the Fermi surfaces are of 
equal sizes. As one increase the temperature the system de¬ 
velops an imbalance in the population of the up and down 
fermionic species, the signature of which is observed in the 
increasing size mismatch between the two Fermi surfaces. 
There is already a weak signature at T = 0.1 It, a clear sig¬ 
nature at T ~ 0.131 ^ T c (not shown here), and a prominent 
difference at T = 0.21 and T = 0.31. 
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FIG. 7. Color online: Thermal evolution of the various indicators at h = 0.5 t. Starting from the top we show maps of | A|, phase correlation, 
pairing structure factor, magnetization m*, magnetic structure factor, and number density. The temperature, along the row, is marked at the 
bottom of the figure. Between T — O.lOt and 0.13t, see Fig.5, there is both significant superfluid order as well as magnetization. 
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FIG. 8. Color online: Thermal evolution of the momentum occupa¬ 
tion number n a ( k) at h — 0.50. The up and down spin distributions 
are same at T = 0, slightly different for T < T c (2nd row), and 
noticeably different at high temperature. 


D. High field: appearance of metastable FFLO states 

In the high field regime, 0.7 1 < h < 0.85 1, the region 
of first order USF to normal transition, the system seems to 
encounter competing minima in the energy landscape. The 
state we obtain depends on the thermal history of the system. 
We highlight the effects at a typical field h = 0.8£. 

Fig.9 shows the standard spatial and Fourier space indica¬ 
tors on a cooling run. The results on heating from the USF 
are qualitatively similar to what we have seen at h = 0.5t. On 
cooling from high T the system encounters 0 fluctuations 
and instead of transiting to a q = (0,0) low T state it actually 
enters a modulated state! This state has higher energy than the 
variational USF state which suggests its metastable character. 
We show the real space patterns at the lowest T further on. 
While real space features are not very illuminating down to 
T = 0.05£, the pairing structure factor shows a clear finite q 
feature. The spatial character becomes clearer at even lower 
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FIG. 9. Color online: Thermal evolution of the various indicators on cooling the system at h = 0.8 1 . From top to bottom we have plotted 
the spatial map of |A|, phase correlation, the pairing structure factor, magnetization, magnetic structure factor and number density. This 
illustrates the emergence of a strong LO signature in the pairing structure factor, when the ground state is actually USF. The pattern at the 
lowest temperature: T = 0.00It is shown later. 


T as we show below. 

We compute the momentum occupation numbers for the up 
and down fermionic species through the heating and cooling 
cycles. Apart from the evolution of the mismatch between the 
up and down distributions with temperature one can also see 
the modification in the Fermi surface shape with respect to 
what one would expect in the simple tight binding case. The 


straight Fermi surface segments at low T in Fig. 10 result from 
line-like LO correlations as we show in the Fig.ll. The rise 
in temperature wipes out this feature. 

What does this metastable LO state look like in real space? 
We computed the amplitude, phase, magnetization and num¬ 
ber density maps for MC snapshots and show a typical set at 
low temperature in Fig. 11. As can be seen, real space periodic 
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FIG. 10. Color online: Momentum occupation numbers n CT (k) at different temperatures through the heating and the cooling cycle computed 
at h = 0.8 1. Notice the ‘size difference’ persisting to low temperature in the cooling run - suggesting a finite m ‘ground state’ (actually 
metastable in this case). For T > T c the heating and cooling results are essentially similar. 



FIG. 11. Color online: Spatial maps characterizing the metastable LO state through (a) |A*|, (b) phase correlation, (c) magnetization and 
(d) number density distribution at T — 0.00It. 


modulations are observed in both the superfluid order param¬ 
eter and local magnetization. The order parameter exhibits a 
nodal, domain wall like structure, in the nodes of which reside 
the unpaired fermions giving rise to a finite magnetization. A 
‘node’ in the |A^| corresponds roughly to a peak in the local 
magnetization. 

Before we end this section we show the ideal momentum 
occupation number n^k) corresponding to the metastable 
state at h = 0.80 1 in Fig. 12. A weaker variant of the same has 
been observed and presented in Fig. 10. Fig. 12 prominently 
shows the anisotropic deformation of the Fermi surfaces in 
presence of an underlying modulated pairing order. 


E. Density of states 

Pseudogap features in the density of states, and momen¬ 
tum resolved spectral functions, have been explored experi¬ 
mentally in the balanced ‘continuum’ Fermi gas at unitarity. 


There is also a body of associated theory 76 78 . Unfortunately 
there are no such detailed spectral experiments in the imbal- 
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FIG. 12. Color online: Momentum occupation function for an ideal 
diagonal stripe phase to mimic the pattern that we observe in Fig. 11. 
Fig.l 1 modulations have a 2D character (rather than simple diagonal 
stripe) so the actual n(k) in Fig. 10 has an approximate fourfold look. 






















FIG. 13. Color online: Temperature dependence of the DOS at h — 0.21,0.5£, 0.8 1 (left to right). The top row shows N^(uj) and the bottom 
row shows the total DOS N(uj). Ni(cj) is just a shifted version of While the main feature in (a) is slow filling up of the gap (with 

increasing T the gap converts to a pseudogap already below T c ), (b) and (c) reveal that at higher fields this ‘filling up’ process is non monotonic. 
We quantify the relevant temperature scales later in Fig. 16. 


anced case. Our coupling corresponds roughly to what would 
be considered ‘unitary’ (see the Discussion section), but we 
are working on a lattice, at densities far from the continuum 
end. To check the usefulness of our approach in capturing the 
qualitative features of this well studied end we compared our 
‘balanced’ results to those in the literature 76 . We found that 



FIG. 14. Color online: (Left) Temperature dependence of the DOS 
at the shifted Fermi level for up spin fermions. Right: a color map 
of the DOS at u — — h for varying h and T. The non monotonic T 
dependence is clearly visible in the large h regime. 


the dispersion and damping share several features (we will put 
this up separately) providing confidence that our lattice results 
would have value in analysing imbalanced continuum gases as 
well. 

Fig. 13 shows the spin resolved and total density of states 
at low, intermediate, and high fields, h = 0.2£, 0.5£, 0.8 t, 
respectively. We focus on the ‘up spin’ DOS, since the ‘down 
spin’ DOS is symmetrically shifted (and the total is simply a 
sum of these two) and define Tow energy’ as cc ~ — h. 

The DOS in the top row in Fig. 13 reveal the following fea¬ 
tures: (i) A hard gap with the usual gap edge coherence peaks 
at low T. (ii) The hard gap converts to a pseudogap at a tem¬ 
perature Tpgl, where Tpgl (h) < T c (h ), and the coherence 
features are suppressed with increasing T. The two features 
above are visible at all fields, see panels (a)-(c). An additional 
feature is visible in (b) and (c). This is (iii) the weakening of 
the PG (i.e, increase in the low energy DOS), continues up to 
a temperature T max (h), beyond which the DOS at low energy 
falls again. The low energy DOS would probably flatten out 
and the pseudogap close at some high temperature. The some¬ 
what exotic look in the total DOS, panels (d)-(f) arises from 
adding up and so we focus on N^(uj) itself. 

Fig. 14 shows the T dependence of the spin up DOS at the 
up spin ‘Fermi level’, uo = —h, in the left panel. The right 
panel shows a map of this density of states as a function of h 
and T. These data allow us to extract the temperature T rnax at 
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FIG. 15. Color online: Field dependence of the N^(lj + h), the total DOS N(u) and P(| A|). Top row: T = 0.05 1, bottom row T = 0.15 1. 
The shift in is to gauge out the field dependent shift of origin and the resulting clutter in the plot. 


which the spin resolved DOS at uj = ±/i has its maximum. 

Fig. 15 shows h dependence at fixed temperatures, high¬ 
lighting the effect on the DOS as the system evolves from the 



FIG. 16. Color online: Temperature scales associated with the be¬ 
havior of the spin resolved DOS. The low T hard gap in the super¬ 
conductor converts to a pseudogap at a temperature T pg \ < T c , while 
the ‘ungapped’ partially polarized Fermi liquid at large h develops a 
pseudogap at T = T pg 2 . The entire window above Tpgl and Tpg2, is 
pseudogapped. The DOS at the center of the pseudogap (cu = d= h) 
shows a maximum at T maa; . 


‘balanced’ situation to the highly magnetized (hence weakly 
paired) state. At T = 0.05£, top row, which is ~ 0.3T°, the 
pseudogap in the DOS vanishes at h ~ t while at T = 0.15t, 
bottom row, the PG seems to persist even at h = 1.5 1 where 
the ground state is an unpaired Fermi liquid! 

In panels (c) and (f) we show P(| A|) for the same field val¬ 
ues as in the DOS panels. At T = 0.05 1 the P(|A|) remains 
almost unchanged for h between ^ 0 — 0.61 The resulting 
DOS also remains essentially unchanged over this field win¬ 
dow. At h = O.St the center of P(|A|) is at a noticeably 
smaller value, and as h increases the peak and the mean value 
of | A | shift to progressively lower value. At h = 1.5t the peak 
value is ~ 0.25A 0 and the P(| A|) cannot generate a pseudo¬ 
gap in the spectrum. 

At T = 0.15£ the peak location of P(|A|) and the mean 
shift to lower value with increasing h. All the data, except at 
h = 0, 0.2£, are at T > T c . However, at this temperature 
the mean value at high fields is significantly larger than what 
we see at T = 0.051 At h = 1.5 1 the peak of P(|A|) is at 
~ 0.5Ao, almost twice the T = 0.05t value. 

As a result, even though the high field system starts at low T 
as essentially an ‘uncorrelated’ partially polarized Fermi liq¬ 
uid (within our scheme) the thermally generated correlation 
effects are strong enough to generate a pseudogap with in¬ 
creasing temperature. The spin resolved DOS at large h starts 
gapless (at T = 0) but transits to an interaction induced weak 
pseudogap phase at high T. This pseudogap has nothing to do 
with long range order in the ground state. 




















































Tracking the field dependence of the pseudogap formation 
scale, starting from the low temperature end, allows us to con¬ 
struct the PG feature based ‘phase diagram’ in Fig. 16. It re¬ 
veals several intriguing features: (i) Although the T = 0 gap 
in the spectrum remains the same for ft, = 0 — 0.85£ the tem¬ 
perature T pg i, at which this gap converts to a pseudogap, col¬ 
lapses as ft ft c i the USF-FFLO boundary, (ii) Although 
the mean field ground state has no pairing for ft > 1.3ft and 
is therefore gapless, fairly modest temperature ^ 0.07 1 gen¬ 
erates an weak pseudogap due to thermal generation of pair¬ 
ing fluctuations, (iii) The PG in the spin resolved DOS sur¬ 
vives to a high temperature, certainly greater than T ~ 0.5£ 
that we have probed, although at large ft it is a weak feature. 
This survival to high T is a consequence of the large interac¬ 
tion U = At that we have chosen, and has a parallel in the 
PG observations made on the imbalanced cold Fermi gas at 

unitaritjEU. 


IV. DISCUSSION 

Till now we have mainly focused on our specific results. In 
what follows we touch briefly on a few broader issues. These 
include: (a) the reliability and limitations of our method, (b) a 
conceptual framework for understanding the numerical data, 
(c) the connection between our intermediate coupling lattice 
results and the unitary continuum gas, (d) qualitative compar¬ 
ison with cold atom and solid state experiments, and (e) the 
wider possibilities of our method in exploring imbalanced su¬ 
perfluids in other situations. 


A. Issues of method 

Our results are based on (i) a Hubbard-Stratonovich decom¬ 
position of the interaction in the pairing channel, (ii) approxi¬ 
mating this auxiliary pairing field A^ as classical, (iii) a cluster 
algorithm based Monte Carlo sampling of the field, (iv) use 
of finite size, as is inevitable in any calculation of this kind, 
(ii), (iii) and (iv) introduce errors and we comment on these in 
the paragraphs below. 


1. Hubbard-Stratonovich decomposition 

The analytic basis of the HS based method is discussed in 
the Appendix. 


2. The static approximation 

The static auxiliary field approximation is exact as T —» oc 
and in principle becomes less and less accurate as T —>> 0 (as 
the energy difference between the bosonic Matsubara frequen¬ 
cies reduce). However , when the ground state has some kind 
of long range order, as in both the balanced and unbalanced 
fermion cases, the static mode succeeds in capturing much of 


the interaction effects. This keeps our T = 0 results qualita¬ 
tively valid. A comparison in the balanced case revealed that 
by the time T ~ T c , the static mode captures most of the ther¬ 
mal effects, and anyway for T T c it should describe the 
problem exactly. Overall, in the current problem, the static 
approximation by itself is not a serious limitation. 


3. Single channel decomposition 

A single field decomposition that is static cannot in gen¬ 
eral capture instabilities in all channels. In the FFLO regime 
the pairing, density, and magnetic channels are in principle 
all relevant. However, we find that for our chosen mean den¬ 
sity, the density modulations in the FFLO phase are very weak 
so ‘density channel’ effects are not important (they would be 
very important if n = 1). The presence of an additional mag¬ 
netic channel may make a quantitative difference to our re¬ 
sults. While these additional channels are readily incorporated 
within MFT a non Gaussian fluctuation theory, like ours, in¬ 
volving all these modes is difficult to construct. We have opted 
to stay with a simple decomposition so that the fluctuation the¬ 
ory can be better handled. 

For the homogeneous BP phase we have found that there 
is no density wave ordering tendency at n = 0.94 and the 
field regime that we have considered. Fluctuations in density 
are present, as evidenced in Figs.7 and 9, but are small since 
the field induced magnetization suppresses the density wave 
susceptibility. In a more elaborate calculation the fluctuations 
could be numerically larger, without changing the qualitative 
features of our result. 


4. Monte Carlo: cluster algorithm and size dependence 

The MC implementation using the Bogoliubov-de Gennes 
scheme requires repeated diagonalisation of the fermion prob¬ 
lem. Done exactly this computation scales as TV 4 where N is 
the system size, limiting one to TV ^ 10x10, hardly ad¬ 
equate to access complex phases. This is a primary limita¬ 
tion in FFLO studies and limits most finite temperature stud¬ 
ies to mean field theory. We can access much larger size (up 
to 40 x 40, say) since we use a cluster based update scheme, 
discussed in the text. Unfortunately the cluster size introduces 
another length scale, that affects access to FFLO phases, but 
does not seem to have much impact on the uniform SC state. 
So, as far as the present study is concerned, size limitations 
have not been significant. We have checked the quality of the 
MC in the ft = 0 problem earlier by comparing to full QMC®. 


B. Landau-Ginzburg framework 

It is useful to put up a Ginzburg-Landau (GL) framework 
for qualitatively understanding our results, focusing on a A^ 
only theory rather than the ‘fermion + Aft problem. At 
weak coupling GL theory could have been systematically 
derived 80 81 , here it serves as a phenomenological construct. 




The free energy density suggested by Casalbuoni et al . 8 ^, 
for the superfluid in the presence of a magnetic field, is: 

7 = ^a|A | 2 + 1/3|A | 4 + |t|A | 6 + e|VA | 2 + ||V 2 A | 2 

The complicated form, involving a 6 th order amplitude term 
and V 2 A, is retained since (3 and e which are positive in the 
ft = 0 case can change sign when ft, ^ 0 . 

In the ft = 0 functional involving only a, (3 and e, we have 
(3 > 0. The sign change of a drives a second order transition 
to a q = ( 0 , 0 ) state since the gradient term penalizes spatial 
modulation. 

(3 changing sign from positive to negative leads to a first 
order transition, again to an uniform state if e > 0, and one 
retains a positive 7. On the other hand if e changes sign the 
system would head towards a modulated state, whose wave 
number has to be decided by the presence of a positive 7 . This 
would be the thermal transition to some FFLO state. 

In the continuum weak coupling limit it turns out that (3 and 
e change sign from positive to negative at the same poinW 
In that situation one has a second order normal to SC transition 
at weak field, crossing over to a first order normal to FFLO 
transition beyond a critical field. 

Our lattice mean field results at U = 4 1 indicate that a first 
order thermal transition need not be necessarily to an FFLO 
state. We do have a window of a first order normal to uniform 
SC transition. This distinction is probably a lattice versus con¬ 
tinuum difference. It shows up in the MC results as well, with 
T c scales suppressed due to amplitude and phase fluctuations. 

The MC results suggest the rough behavior of the various 
GL coefficients in terms of ft and T but the observed metasta¬ 
bility is harder to pin down. For that a more elaborate func¬ 
tional, derived by tracing out the fermions from the coupled 
problem, and expanded about q = 0 and q = Q (the FFLO 
wave vector) would be needed. For ft < h\ that free energy 
has the deepest minimum at q = 0 , and also uniquely reaches 
this state on thermal cycling. For fti < ft < h c \ the absolute 
minimum is still at q = 0 but it seems that the minimum at Q, 
although metastable, dominates the energy landscape. When 
one cools from the high T state the system seems to first en¬ 
counter this q / 0 minimum and tracks this state down to 
T = 0. 


C. Connection to continuum unitary gas 

While we have motivated our lattice model in terms of ex¬ 
periments on the continuum unitary Fermi gas there are issues 
which need highlighting. These are (i) the notion of ‘uni- 
tarity’ in the 2D lattice model (and its relation to BCS-BEC 
crossover), and (ii) the access to continuum ‘universal’ effects 
via lattice simulations. 

Unitarity: The primary scale quantifying the strength of 
interaction in cold Fermi gases is the two body 5 -wave scat¬ 
tering length: ap>, where D denotes the spatial dimensional¬ 
ity. The dimensionless coupling constant can then be written 
in terms of kpciD, where k F is the Fermi wavevector. We 


quickly comment on unitarity and the BCS-BEC crossover is¬ 
sue in the continuum and lattice contexts and then see what 
information lattice simulations can yield. 

In the 3D continuum Fermi gas the scattering length a^D 
oc as the interaction g —>> g c . g c is a finite in 3D, and 
defines the interaction for which a two body bound state 
first forms in vacuum. The (inverse) dimensionless coupling 
1/kpasD = 0 at g c . This is also the point near which the tran¬ 
sition temperature of the 3D Fermi gas has its maximum, with 
T™ ax /E F ^ 075 On the 3D Hubbard lattice an equivalent 
critical interaction for bound state formation can be worked 
out and yields U c /t ~ 7.9. Again, the T c is found to be maxi¬ 
mum for U/t ~ 8 (from QMC), remarkably closed to U c . 

In the 2D continuum a two body bound state forms in the 
presence of an arbitrary attractive interaction so < 221 ? (g) —> oc 
for g 0. This is however in the deep BCS regime where 
a weak coupling description in terms of fermionic quasiparti¬ 
cles is sufficient. As g increases the pair size shrinks and in the 
Bose limit a^D 0. The crossover coupling in the 2D case 
is defined via ln{k F a 2 D) 0 , i.e, the scattering length being 
comparable to inter particle separation. We are not aware of a 
2D continuum QMC calculation for the T c , but interpolation 
between the BCS and BEC end suggests that the maximum T c 
occurs for ln(k F d 2 D ) —>* 0, wit1oMT™ ax /E F ~ 0.1. 

On the 2D Hubbard lattice a two body bound state would 
form at arbitrary weak attraction, i.e, U c /t —» 0. The notion 
of fti? is not very meaningful on the lattice, particularly away 
from low density, but following the cases above one may iden¬ 
tify the crossover as the region of maximum T c . The T c is well 
established via QMC and the maximum occurs for U/t ~ 5. 
For U/t = 4 that we use the T c is ^ 0.9 of the maximum 
value 73 . 

Overall, the two body bound state based unitary point in 3D 
also corresponds to a regime where (i) neither the fermionic 
nor the bosonic quasiparticle description suffices, and (ii) the 
T C /E F is maximum. In 2D the simple two body argument 
would put this regime at extreme weak coupling but (i) and (ii) 
above indicate that [ln(k F a 2 D)]~ 1 oc, rather than a^D 
00, is the relevant choice. 

Our interaction strength is not far from what would be con¬ 
sidered the ‘unitary’ value in the 2D lattice case. Can we 
comment on the universal physics one would have seen in the 
continuum case, on which, remarkably, there is now an exper¬ 
iment? 

Continuum universality from lattice physics: The first 
difficulty is with our density choice, we have used n ^ 0.9 
to maximize T c (at the same time avoiding the density wave 
instability at n = 1). At this density the lattice effects are 
very prominent and the Fermi surface is distinctly non cir¬ 
cular. Even if we had used lower density, n ~ 0.1, say, 
where the Fermi surface is indeed circular and ek ~ ft 2 is 
a good approximation access to continuum effects is difficult. 
If the interaction were weak, i.e, U <C 4t, the physics would 
have been insensitive to the high energy band cutoff. How¬ 
ever, ‘unitarity’ requires U ~ 5t and even if the Fermi level 
is at the lower edge of the band, scattering effects couple in 
states at the upper edge. The high energy states are lattice 
specific, and as the paper by Privitera et aP * 85 demonstrates 





FIG. 17. Color online: The imbalance-temperature phase diagram 
inferred from measurements on a cold atomic gas at unitarity (leftj^j, 
compared to our result on the intermediate coupling (peak T c ) Hub¬ 
bard model (right). The normalization of the x axis is same in both 
panels, while the y axis have different reference scales. 


even at reasonably low densities one obtain results that do not 
match with continuum predictions. So, extracting ‘universal’ 
physics from lattice simulations require extremely low densi¬ 
ties, ^ 0 . 001 , or lattice sizes ~ 10 4 , outside the range of what 
is doable today. 


D. Comparison with experiments 

1. Atomic superfluids 

Our model finds it most appropriate experimental coun¬ 
terpart in imbalanced cold Fermi gases, studied at strong 
attractive interaction promoting s-wave pairing. There are 
still differences, e.g , (i) our results are on a lattice theory, 
the experiments are in the ‘continuum’, (ii) there is a trap 
present in the experiments, and (iii) dimensionality (the ex¬ 
periments are in three dimensions). Nevertheless, the simi¬ 
larities are striking. Fig. 17 presents the experimental phase 
diagram 32 of the unitary Fermi gas in terms of magnetization 
m = (n^ — n±)/(n^ + n±) and temperature, constructed by 
the experimental group by ‘gauging out’ the effect of the trap. 
Next to it we show our m — T phase diagram, constructed at 
n ~ 0.9 by varying h (hence m). 

The experiments infer (i) a homogeneous magnetized su¬ 
perfluid (SF), (ii) an unstable region, and (iii) a magnetized 
normal Fermi liquid. At T = 0 the SF to unstable transition 
occurs at m = 0+, suggesting that a finite m homogeneous 
SF state cannot occur at T = 0, while the ‘unstable’ to normal 
transition occurs at m ~ 0.35. The T c of the unitary Fermi gas 
is significantly lowered with respect to mean field theory 82 !and 
using the measured zero imbalance T c as the reference scale, 
the tricritical point occurs at m ~ 0.2 and T t ri ~ 0.4T C °. The 
suppression of T c in presence of imbalance is an extension 
of the zero field (balanced) behavior 51 53 55 . The presence of 
imbalance makes the suppression rapid. 

The m — T picture that emerges from our data already has 


the fluctuation effects built in on T ®. In the ground state the 
SF to ‘unstable’ transition occurs at m = 0+ and the unstable 
to LO transition at m ~ 0.28 and the LO to normal transition 
at m ~ 0.37. It must be noted that the unstable region es¬ 
sentially is a phase separated region, marked by discontinuity 
in density, and characterized by the absence of homogeneous 
superfluid phase. The LO state has not been observed ex¬ 
perimentally in the 2D geometry, possibly due to additional 
fluctuations in the absence of a lattice. Our tricritical point is 
at m ~ 0.15 and T tri ~ 0.07t ~ 0.4T,?. Given that we are in 
a regime where the Fermi surface is significantly non circular 
the overall correspondence of phase boundaries and tempera¬ 
ture scales is reasonable. We have predictions about spectral 
properties on the m — T plane that we will present separately. 


2. Superconductors 

The solid state systems in which Pauli limited behavior 
is observed, for example CeCoIn 5 and the organic super¬ 
conductors, are non s-wave materials and have Tow energy’ 
fermionic degrees of freedom even in the ordered state. In our 
model, however, the fermions are gapped, or have a strong 
pseudogap, due to the large on site attraction - unless the pop¬ 
ulation imbalance is large. As a result, despite the overall 
similarity in the look of our theory phase diagram and those 
observed in experiments®^, the comparison of indicators 
like specific heat, Cy (T, h), and magnetization, m(T, h), re¬ 
veal differences in detailed behavior. We have made these 
comparisons but do not present the data here. 

E. Extensions of the present method 

The present work was focused on understanding a part of a 
larger phase diagram. As a natural extension of this we have 
studied the thermal properties of the large h FFLO states in de¬ 
tail. We have also computed the momentum resolved spectral 
functions of the BP, PPFL, and FFLO phases over the entire 
h — T window. We will present these results separately. 

A natural extension of the present method, involving a ‘two 
field’ decomposition, can handle the effect of disorder^® on 
the FFLO state, including the thermal effects which are in gen¬ 
eral difficult to access. 

Finally, cold Fermi gases involve a trapping potential and a 
non trivial spatial dependence of the region where the fluid 
is magnetized. While experimental optical lattice sizes ^ 
100 x 100 are hard to access using our MC technique, we 
hope to access the physics at least in the BP regime using a 
local density scheme grafted on to our Monte Carlo solver. 


V. CONCLUSION 

We have used a real space Monte Carlo technique based on 
a static pairing field approximation to study the behavior of 
a Pauli limited superconductor in the BCS to BEC crossover 
regime. We find that the T c scales are strongly suppressed 








with respect to mean field predictions, there is a wide window 
of metastable FFLO states in which the system gets trapped 
when the true ground state is a homogeneous superfluid, and 
the spin resolved density of states shows a non monotonic low 
energy character. We do not know of Pauli limited solid state 
systems with s-wave pairing, but ultracold unitary gases sug¬ 
gest an universal phase diagram quite similar to what we ob¬ 
serve. This paper probes the lower field ‘breached pair’ state 
in detail, companion papers discuss the FFLO regime and the 
spectral features expected with changing imbalance. 
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VI. APPENDIX: HUBBARD-STRATONOVICH 
TRANSFORMATION AND MONTE CARLO SAMPLING 


The primary numerical technique we use is a Monte Carlo 
implementation of a ‘single channel’ static auxiliary field de¬ 
composition of the A2DH]V ^ 9 I — I l 89 l 9Q l Below we discuss about 
the various aspects of our numerical technique. 

The Hubbard model at strong interaction requires a non per¬ 
turbative solution. The exponential growth in the dimension 
of the Hilbert space rules out the use of exact diagonalization 
except for very small sizes. The ‘exact’ tool of choice is quan¬ 
tum Monte Carlo (QMC) against which all approximations are 
bench marked. While QMC can be implemented via various 
approaches, the method below easiest reveals the connection 
to our approach. 

The Hubbard partition function is written as a functional 
integral over Grassmann fields ^ CT (r), 

Z = J VipV^e-^’^ 

rP 

S= dT\S2{ip ia (d T 5ij + tij)ip ja } - \U\ V V’ifV’itV’uV’tl.] 

ij,<T i 

Only quadratic path integrals can be exactly evaluated. Since 
the interaction generates a quartic term in the - 0 ’s the partition 
function cannot be immediately evaluated. 

The quartic term is ‘decoupled’ exactly through a Hubbard- 
Stratonovich transformation in terms of pairing fields 
Ai(r), Ai(r). This introduces a term A^^(t)' 04 (t) in the 
action. 

Z = jv AVA*V^V^e~ Sl[ ^ AA * ] 

rP 

S\ I dT^^^\jpi a {d T 8ij ~b ) 'ipj a } 

° ij,c r 

+Lb Ai ( T )^t( T )^4-(j)+ h - c +^j-i] 


The integral is now quadratic but an additional integration 
over the field A fir) has been introduced. The ‘weight factor’ 
for the A configurations can be determined by integrating out 
the ip, Vb an d using these weighted configurations one goes 
back and computes fermionic properties. Formally 


■/ 


VAVA*e~ S2[AA * ] 


S 2 = log[Det[g- L - A]] + 


\U\ 


where Q is the Greens function associated with the non inter¬ 
acting H. 

The weight factor for an arbitrary space-time configuration 
A i(r) involves computation of the fermionic determinant in 
that background. If we write the auxiliary field A$ (r) in terms 
of its Matsubara modes, as A*(fi n ), then the various approxi¬ 
mations can be readily recognized and compared. 


• Quantum Monte Carlo retains the full ‘z,fl n ’ depen¬ 
dence of A computing log[Det\^~ x — A]] iteratively 
for importance sampling. The approach is valid at all 
T, but does not readily yield real frequency spectra. 

• Mean field theory restricts A i(Q n ) to a spatially uni¬ 
form (or periodic) and time independent (Q n = 0) 
mode, i.e, A i(iCl n ) —» A. The free energy is mini¬ 
mized with respect A. When the MF order parameter 
vanishes at high temperature the theory trivializes. 

• Our static auxiliary field (S AF) approach retains the full 
spatial dependence in A but keeps only the D n = 0 
mode, i.e , A i(£l n ) —>> A im It thus includes classi¬ 
cal fluctuations of arbitrary magnitude but no quantum 
(Q n 7 ^ 0) fluctuations. One may consider different 
temperature regimes: (1) T = 0: since classical fluc¬ 
tuations die off at T = 0, SAF reduces to standard 
Bogoliubov-de Gennes (BdG) MFT. (2) At T / Owe 
consider not just the saddle point configuration but all 
configurations following the weight e~ S2 above. These 
involve the classical amplitude and phase fluctuations of 
the order parameter, and the BdG equations are solved 
in all these configurations to compute the thermally av¬ 
eraged properties. This approach suppresses the or¬ 
der much quicker than in MFT. (3) High T: since the 

= 0 mode dominates the exact partition function 
the SAF approach becomes exact as T oo. 


• DMFT: for completeness we mention that DMFT re¬ 
tains the full dynamics but keeps A at effectively one 
site, i.e , A^(£ 2 n ) y A (£ 2 n ). 


Overall, our method reduces to BdG mean field theory only 
at T = 0 but retains all the classical thermal fluctuations at 
T / 0. As a result it is only as good as MFT at T = 0 
but is far superior in estimating T c , and essentially exact as 
T —y oo. It does use BdG iteratively as a tool but on all 
fluctuating configurations not just the mean field state. 
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